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Abstract 


In 1999, Stolz and Adams [Phys. Fluids 11 (1999)] unveiled a subgrid-scale model for LES 
based upon approximately inverting (defiltering) the spatial grid-filter operator and termed 
.the approximate deconvolution model (ADM). Subsequently, the utility and accuracy of the 
ADM were demonstrated in a posteriori analyses of flows as diverse as incompressible plane- 
channel flow and supersonic compression-ramp flow. In a prelude to the current paper, a 
parameterized temporal ADM (TADM) was developed and demonstrated in both a priori 
and a posteriori analyses for forced, viscous Burger’s flow. The development of a time- 
filtered variant of the ADM was motjvated*primarily by the desire for a unifying theoretical 
and computational context to encompass direct numerical simulation (DNS), large-eddy 
simulation (LES), and Reynolds averaged Navier-Stokes simulation (RANS). The resultant 
methodology was termed temporal LES (TLES). To permit exploration of the parameter 
space, however, previous analyses of the TADM were restricted to Burger’s flow, and it has 
remained to demonstrate the TADM and TLES methodology for three-dimensional flow. 
For several reasons, plane-channel flow presents an ideal test case for the TADM. Among 
these reasons, channel flow is anisotropic, yet' it lends itself to highly efficient and accurate 
spectral numerical methods. Moreover, channel-flow has been investigated extensively by 
DNS, and a highly accurate data base of Moserebal. [Phys. Fluids 11 (1999)] exists. In the 
present paper, we develop a fully anisotropic TADM model and demonstrate its utility in 
simulating incompressible plane-channel flow at nominal values of Re r = 180 and Re r = 590 
by the TLES method. The TADM model is shown to perform nearly as well as the ADM at 
equivalent resolution, thereby establishing TLES as a viable alternative to LES. Moreover, 
as the current model is sub-optimal is some respects, there is considerable room to improve 
TLES. 
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I. INTRODUCTION 


While the formal linkage of the equations governing large-eddy simulation (LES) and 
Reynolds-averaged Navier-Stokes simulation (RANS) has been well established, 9,10 it is of 
interest to investigate whether this linkage can be extended practically by developing fil- 
tering and averaging procedures that yield mutually consistent solution fields among DNS, 
LES, and RANS. A possible unifying context for these methodologies is afforded by filter 
theory. However, the linkage between DNS, LES, and RANS may be more natural within 
the context of time-domain filtering rather than the traditional spatial filtering commonly 
used in LES. Accordingly, the present paper attempts to establish temporal large-eddy sim- 
ulation (TLES) as a practical methodology for solving the temporally filtered Navier-Stokes 
(TFNS) equations using causal time-domain filters. TLES lies at the nexus of three rela- 
tively recent developments in subgrid-scale (SGS) modeling for LES: dynamical modeling, 
approximate deconvolution methods, and time-domain filtering. Each is discussed in some 
detail below. 

A. Dynamic Modeling 

LES dates to the early 1960’s, when it was first exploited for weather modeling. Summa- 
rizing Pope 22 , LES consists of four conceptual steps: 1) decomposition of the field variables 
into resolved and unresolved scales of motion, denoted here by u(t,x ) and u(t, x), respec- 
tively; 2) derivation of the equations of motion for the resolved scales; 3) closure of the 
governing equations by modeling the residual-stress tensor; and 4) numerical solution of the 
closed governing equations. Only in the last step does discretization arise. 

Steps 1) and 2) are accomplished by subjecting the flow-field variables and the Navier- 
Stokes Equations (NSE), respectively, to the same low-pass filter, conventionally termed the 
grid filter. The terminology is unfortunate, because the filter is applied prior to discretiza- 
tion; hence, we will refer instead to the primary filter. Filtering the nonlinear terms of 
the NSE generates residual stresses, which require modeling in step 3. The more common 
term for residual stress is subgrid-scale (SGS) stress, but for the reason stated previously, 
this too is misleading. More correctly, the residual stress is alternately referred to as the 
subfilter-scale (SFS) stress, a distinction made by Gullbrand 13 . In conventional LES, the 
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primary filter is implicit; that is, it is conceptual and does not play an active role in the 
solution procedure (Step 4). 

The original residual-stress models were of eddy-diffusivity type, of which the Smagorin- 
sky model is the most well known. Eddy-diffusivity models suffer from a number of defects, 
among them the tendencies to be overly dissipative and to correlate poorly with exact resid- 
ual stresses. A more subtle shortcoming of initial approaches to LES was that implicit 
filtering obscured the relationship between the filter and the model, so much so that, unt il 
relatively recently, it was widely accepted within the LES community that the choice of filter 
and the model were completely independent 21 ’ 26 . 

In an effort to overcome the shortcomings of conventional LES, Germano and coworkers 8 
proposed the concept of dynamic modeling, which re- invigorated the LES community. Dy- 
namic modeling involves two filters: the primary (“grid”) filter, which is implicit, and the 
secondary ( “test” ) filter, which is explicit. The Germano identity links the residual-stress 
tensor with two other tensors, one of which the resolved turbulent stress tensor, is 
computable by secondary filtering of the resolved velocity fields. The quantity C tJ can be 
considered a measure of ill-resolution, from which a localized dissipation parameter is de- 
rived. Thus, dynamic models axe local in time and space in that dissipation is applied only 
when and where needed. Perhaps the greatest contributions of dynamic modeling are ex- 
plicit filtering and local dissipation; however, dynamic models suffer some theoretical and 
practical shortcomings. 

The advent of explicit filtering resulted in a careful examination of the relation- 
ship between the filter(s) and the model from several points of view: experimental 19 , 
computational 24 , and analytical 25 ’ 30 . All recent investigations concur that there is strong 
coupling between the filter and the exact residual stress. Fax from being independent, the 
model and filter axe intricately related to the point that there is one-to-one correspondence 
between the filter and the model 24,30 . A secondary motivation for the current work is to 
render that dependence explicit. 

B. Approximate Deconvolution Methods 

Conventional LES has relied upon models that are phenomenologically based; that is, 
they axe derived primarily on the basis of physical considerations. Recent recognition of the 
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tight coupling between filter and model has spurred attempts to approximate the residual 
stress mathematically rather than to model it. Taylor-series analyses (Leonard 18 , Bardina 4 , 
Pruett et al 25 ) fall into this category, as do more recent deconvolution methods (Shah and 
Ferziger 29 , Geurts 11 , Domaradzki and Saiki 7 , and Stolz and Adams 31 ). A lucid review of 
deconvolution methods can be found in Adams and Stolz 2 . 

Deconvolution methods are based upon deconvolving (defiltering) the filtered flow fields 
(provided that the filter is invertible). Interestingly, the method enjoys its utility because 
it is approximate, not despite being approximate. In LES, one wishes to recover accurately 
only the resolved scales, which deconvolution enhances, while appropriately dissipating the 
energy of the unresolved scales, which the ADM method also accomplishes by secondary reg- 
ularization. The (spatial) approximate deconvolution model (ADM) of Stolz and Adams 31 
has performed remarkably well in o posteriori analyses of flows as diverse as incompressible 
plane-channel flow 32 and supersonic compression-ramp flow 33 . Moreover, the dissipation 
provided by the filter renders the method applicable for shock-capturing in high-speed com- 
pressible flows (Adams 1 ). 

C. Time-Domain Filtering 

As mentioned previously, LES has relied historically on spatial filtering to separate re- 
solved from unresolved scales. However, at least in concept, time-domain filtering offers 
a number of advantages, which were outlined previously in Dakhoul and Bedford 6 and in 
Pruett 23 . For completeness, we compile some of the advantages below. 

1. Time-domain filters naturally commute with differentiation operators; commutation 
error, however, is problematic for spatial filters, particularly on finite domains 5,35 . 

2. Spatial filtering is problematic for highly stretched meshes or unstructured grids; tem- 
poral filters, on the other hand, operate independently of spatial discretizations. 

3. LES based upon temporal filtering (=TLES) permits consistent comparison with re- 
sults from most physical (e.g., wind-tunnel) experiments, in which data are typically 
recorded and processed in the time domain. 

4. Time domain filters are compatible with flow manipulation by localized (point) 
sources 6 . 
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5. As will be shown, the reformulation of LES for time-domain filters (=TLES), yields 
a parameterized system of governing equations for which the temporal filter- width A 
appears as an explicit parameter. 

6. Last and perhaps most important, time-domain filtering provides a natural and uni- 
fying context to encompass DNS, LES, and RANS methodologies. 

Time-domain filtering for LES also presents some conceptual and practical disadvantages as 
well, which axe deemed relatively minor and are discussed in context. 

Surprisingly little effort has been expended toward the exploitation of time-domain filter- 
ing for LES. Dakhoul and Bedford 6 and Aldama 3 each considered mixed space-time filtering 
for LES. In recent years, attention has been given to pure time-domain filtering, and both 
Eulerian and Lagrangian temporal filters have been considered in lieu of spatial filters for 
LES. In 1996, Pruett 23 applied Eulerian time-domain filtering for a priori analysis of axisym- 
metric jet flow, and in 1999, Meneveau et al. 20 exploited a Lagrangian time-domain filter 
for LES and demonstrated the time-filtered method for LES of isotropic turbulence. In a 
prelude to the current paper 27 , a parameterized temporal ADM (TADM) was developed and 
demonstrated in both a priori and a posteriori analyses for forced, viscous Burger’s flow. 
The analyses were necessarily restricted to Burger’s flow so that the parameter space could 
be fully explored. However, for the spatially one-dimensional problem, LES results agreed 
very well with temporally filtered DNS results, suggesting that TADM methodology should 
be further pursued, which is the intent of this paper. 

In the present paper, we develop an anisotropic TADM model and demonstrate the model 
in large-eddy simulation of incompressible plane channel flow at nominal Re r = 180 and 
590. In the next section, causal filtering is discussed, a candidate filter is presented, and a 
differential form of the filter is derived. In Section III, the temporally filtered Navier-Stokes 
equations are presented. Section IV is devoted to nomenclature. Deconvolution methods are 
discussed in general in Section V, and a TADM model is developed. Results of a reference 
DNS of channel flow are presented in Section VI. Section VII presents results of TLES with 
the TADM and compares these results to results of the reference DNS. Relevant discussion 
follows in Section VIII. The paper closes with few concluding remarks and suggestions for 
future study in Section IX. 
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II. CAUSAL FILTERING 


Time-domain filters are classified as causal or acausal , 34 depending upon whether they 
are applicable to real-time or a posteriori data processing, respectively. The interest here 
lies in real-time applications to TLES for which only causal filtering is appropriate. 

Let /(f) be a continuous function of time t. A causal linear filter is readily constructed 
by the integral operator 

/(*; A )= [ <?(r - t ; A )/(r)dr, (1) 

J -DO 

where G is a parameterized filter kernel , and the parameter A is the termed the filter width. 
(The convention of using semicolons to separate parameters from independent variables in 
argument fists is adopted here.) 

The following properties of admissible kernels were discussed in 27 , but are repeated here 
for completeness. 

l /A 

( 2 ) 




where g is any integrable function such that 

rO 


g(t ) > 0, f g{t)dt = 1 and p(0) = 1. 

J-0 O 

The non-negativity and normalization constraints in Eq. (3) imply that 


fim git) = 0, 

t->-Z > O 


( 3 ) 


( 4 ) 


and suffice for G to approach a Dirac delta function as its parameter A — + 0; that is, 

fim/(f;A) = fim J G(t - f; A)/(r)dr 


6(t — t)/(r)dr 


= m. 


( 5 ) 


In 27 , two examples of causal filters satisfying the constraints above were presented: a 
Heaviside filter and an exponential filter. Here, we restrict attention to the exponential 
filter (for reasons soon to be demonstrated). For the exponential kernel 

exp (t/A) 


g(t) = exp(t) G{t; A) = 


( 6 ) 


( 7 ) 


and the resulting integral operator in Eq. (1) is 

fit A ) = X / 6XP (~A~) f^ dT - 

A drawback of the integral formulation just presented is the need to retain the long-time 
history of the solution field. However, by considering instead the differential form of the filter 
operator, storage requirements are reduced significantly, subject to the intrinsic storage needs 
of the numerical time- advancement scheme itself (for example, low-storage Runge-Kutta). 
By using Leibniz’ rule to differentiate Eq. (7) with respect to time, the differential form of 
the exponential filter is obtained, namely 


d_ 

dt 


/(i;A) = 


/(*)-/(*; A) 


( 8 ) 


The effect of a causal filter is most apparent from its transfer function If (£2), which 
quantifies its amplitude and phase effects in Fourier space as a function of dimensionless 
frequency 0 = u>A, where u is the (dimensional) circular frequency. The transfer function 
of the continuous exponential filter is shown in Fig. 1. When causal filtering is applied to 
a temporally discretized problem with a time increment of At, the action of the filter is 
naturally parameterized by the filter-width ratio r defined as 



(9) 


For the exponential filter, the parameterized transfer function is 

H(vAt;r) = — — (10a) 

(10b) 

Figure 2 presents the modulus of the transfer function of the exponential filter for selected 
values of the filter-width ratio. Note that wA t — tt corresponds to a sampling rate at the 
Nyquist frequency, and that filtering at Q > tt is disallowed because it results in unacceptable 
aliasing error. Note also that r = 0 yields H(uAt\ 0 ) = 1 , which eliminates the filter. 

The filter-width ratio, r, is the only parameter of the differential filter. In general, 
the larger the value of r, the more dissipative the filter. (In this context, a “dissipative” 
low-pass filter is one with significant and broad-band attenuation of high-frequency Fourier 
harmonics.) The differential equation remains viable for all values of filter-width ratio (0 < 
r). However, whenever rs;0, the evolution equation becomes stiff, and small time steps are 
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necessary for stability of the numerical integration scheme. The action of the parameterized 
differential filter on a spectrally rich time signal can be found in Fig. 5 of Reference 27 . 

Normally, when a continuous low-pass filter is discretized, the transfer function of the 
discrete filter deviates substantially from that of the continuous prototype, particularly at 
high frequencies. This is not the case, however, when the continuous filter is expressed in 
differential form, and the differential form is integrated in time by a high-order numerical 
scheme. For example, Fig. 3 compares the transfer function of the continuous exponential 
filter for r = 1 with that of the discrete filter obtained by solving Eq. 8 by the classical 
fourth-order Runge-Kutta method with a time-step At that exactly satisfies the Nyquist 
criterion at the highest frequency. The continuous and discrete transfer functions are nearly 
indistinguishable. Surprisingly, the agreement improves as the filter-width ratio r increases. 
Thus, for present purposes, it suffices to base subsequent analyses on the properties of the 
continuous differential filter. 


III. TEMPORALLY FILTERED NAVIER-STOKES EQUATIONS 


The application of a causal temporal filter (e.g., Eq. (1)) to the Navier-Stokes equations 
(NSE) leads to the temporally filtered Navier-Stokes Equations (TFNSE), namely 

5t-o 

OXj 

dUi t d(uiUj ) _ dp 1 d?Ui dRij 
dxj dxi Redxjdxj dxj ’ 


+ 


( 11 ) 

(12) 


dt dx, dx. 

where Uj is the velocity, p is the pressure, and Re is the Reynolds number. An overbar 
denotes a quantity that has been subjected to filtering by the primary filter, which is causal, 
and R{j represents the exact temporal residual-stress tensor defined as 


Rij — ILitlj IL-iUj . 


(13) 


Provided that filtering and differentiation operations commute, the TFNS equations are 
formally identical to the spatially filtered Navier-Stokes equations. As pointed out previously 
by Pruett, 23 commutativity is natural for temporal filters but remains problematic for spatial 
ones 0,35 . In general, for spatial or temporal primary filters, the residual stress depends 
strongly upon the filter, particularly upon its filter width and order property, which influence 
both the magnitude and the distribution of the residual stress. Because the exact residual 
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stress depends upon the filter width, it is sometimes helpful to explicitly denote by Rj = 
Rij( A). In particular, it can be shown by Taylor-series analysis (e.g., 2a ), that Rij is of 
leading order A 2 for first- or second-order primary filters. 

The invariance properties of the TFNSE are discussed in detail in Section Ilia of 
Reference 27 . 

IV. NOMENCLATURE 

Because TLES represents a departure from the norm, it is advisable to define all terms 
and quantities precisely in the current context, which is the purpose of this section. To those 
for whom the discussion may seem pedantic, we apologize. 

A. Ensemble Mean and Long-Time Average 

Let E{ui} denote the ensemble mean (or the expected, value ) of the velocity component 
U{. Accordingly, the Reynolds decomposition of the velocity is given by 

Ui = E{ui) + u'i ( 14 ) 

which partitions the velocity into time- independent and time-dependent contributions de- 
noted the mean and fluctuation, respectively. Turbulence modeling efforts focus on the 
Reynolds-stress tensor = E{u'u'-}, which, by virtue of the Reynolds decomposition, is 
given by 

Tij = E{(ui - E{ui}){uj - E{uj})} ( 15 ) 

By the ergodic hypothesis, for a statistically steady (stationary) flow, the ensemble mean 
of a turbulent quantity is equivalent to its long-time average, denoted hereafter by angle 
brackets. For stationary flow, for example, < Ui >= E{ui}. In principle, the long-time 
average of a quantity is a constant computed by averaging over infinite time. In practice, 
it suffices to average over many integral time scales (t). If the flow is homogeneous, the 
temporal interval necessary to compute the long-time average is reduced significantly by 
averaging also in homogeneous dimensions. As will be shown, channel flow is homogeneous 
in streamwise (x) and spanwise (y) extents. Consequently, in present, parlance, <> denotes 
either <> t or <>t, x ,y, where the time interval in either case is chosen sufficiently long to 
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yield an accurate mean. In contrast, <> x , y denotes instantaneous averages over homogenous 
planes. 

For channel flow, all field variables are either statistically symmetric or antisymmetric 
about the mid-plane. Due to the random influences of initial data, simulations produce 
slight asymmetries. However, because upper and lower channel halves represent different 
realizations, it is appropriate to average symmetrically about the channel mid- plane. Sym- 
metric averages further reduce the time interval necessary to acquire an accurate mean and 
are exploited herein. 

Henceforth, we presume that the flow is fully turbulent and stationary, in which case the 
Reynolds stress is defined alternately as 

Tij = < ( Uf- < Ui > ) (Uj — < Uj >) > (16) 

= < TJLiUj > — < Ui >< Uj > (17) 

■where Eq. 17 exploits the properties that << a >>=< a >, < ca >= c < a >, and 
< a + b >=< a > + < b > for any time-dependent fields a and b and time-independent field 
c. 


B. Resolved and Sub-Filter Scales 

For TLES, a natural decomposition of the velocity field u t is given by 

Ui = Ui — Ui (18) 

where the first and second terms on the right-hand side of Eq. 18 are termed the resolved and 
sub-filter scale (SFS) velocities, respectively. (Because the TFNSE are continuous and have 
yet to be discretized, we prefer the terminology sub- filter scale of Gullbrand 13 rather than 
the conventional but somewhat misleading term sub- grid scale (SGS).) It is worth noting, 
that, in general, unlike < tij >, Uj is time-dependent. 

In order to formally link RANS and TLES methodologies, we now wish to examine 
the relationship between the Reynolds stress (Eq. 17) and the (temporal) residual stress 
(Eq. 13). The reader is reminded that, unlike the Reynolds stress, the residual stress is a 
time-dependent function, parameterized by the filter width A. 
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In TLES the solution variables are the resolved velocities 5*. Consequently, the exact 
Reynolds stress is unavailable from the solution; the analogous computable quantity is the 
resolved- scale Reynolds stress fy defined as follows: 

fy = < (itj- < Ut >)(Uj- < Uj >) > (19) 

= < UiUj > — < Uj >< Uj > (20) 


C. Exact TLES 

As a sort of thought experiment, we now consider exact TLES; that is, TLES conducted 

with an exact residual-stress model. In exact TLES, the Reynolds stress is exactly the sum 

of the resolved-scale Reynolds stress and the mean residual stress; i.e., 

Ty = < UiUj > - < Ui >< Uj > 

= < TRUj > - < Ui >< Uj > ( 21 ) 

= [< TRUj — UiUj >] + [< UiUj > — < Ui >< Uj >] 

= < Rij > +Ty 

Of the two contributions to the exact Reynolds stress, denoted above by bracketed terms, 
the first, the resolved-scale Reynolds stress, is computable from the solution, whereas, the 
second, the residual stress, must ultimately be modeled. For the moment, we presume the 
residual- stress model to be exact. Still, it may not be immediately clear why Eq. 21 is 
identical to Eq. 17. First, it is useful to note that Eq. 21 is exact in the limiting cases as 
A — + 0 and A — ► oo. By virtue of Eq. 5, Uj — » Ui and HRZJ — * UiUj as A — ► 0. Hence, in 
the limit of vanishing A, the mean residual stress (bracketed term 2 of the second line of 
Eq. 21) vanishes as the resolved-scale Reynolds stress (term 1 above) tends toward the exact 
Reynolds stress. On the other hand, it was shown in Reference 27 that the residual stress 
tends asymptotically toward the Reynolds stress as A — ► oo; that is, limA -,00 Rij = r ij for a 
stationary flow, this by virtue of the fact that limA-.oo /(t; A) =< /(£) > for causal low-pass 
filters. Consequently, in the infinite limit, the first bracketed term in Eq. 21 tends toward 
« Ui >< Uj >> — < u^ >< Uj >= 0; that is, the resolved-scale Reynolds stress vanishes. 
For finite A, the last line of Eq. 21 is also identical to Eq. 17, for the reason that, in exact 
TLES, < ui >=< Ui >, because low-pass temporal filtering preserves the mean. 
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For conventional (spatial) LES, it is customary to validate the results against filtered DNS 
results; that is. by a posteriori analysis. For spatial LES. it suffices to spatially filter the 
DNS results at instants in time and then average over relatively few time steps. In contrast, 
a posteriori analysis for TLES requires real-time or ex post facto temporal filtering of DNS 
data over many time steps. Thus, to conduct a posteriori analyses for TLES, one would 
be forced into DNS with the substantial computational overhead associated with real-time 
implementation of the SFS model or with the enormous storage overhead associated with 
storing most or all time steps. Fortunately (provided the SFS model is adequate), Eq. 21 
obviates the need for a posteriori analysis altogether, because it provides a mechanism by 
which to compare TLES results directly with DNS results without the need to filter the 
latter. Indeed, this is an advantage of TLES relative to spatial LES, because conventional 
LES provides no means for directly comparing LES and DNS or LES and experiment. 

To summarize, the situation is thus. In TLES, the Reynolds stress is partitioned into 
two parts: the resolved-scale Reynolds stress f, ; , which can be computed from the solution 
of the TFNSE, and the residual stress Rij, which must be modeled. The exact partitioning 
depends on the temporal filter width A. For vanishing filter width, all the Reynolds stress 
is contained in resolved scales, the residual stress vanishes identically, the TFNSE reduce to 
the NSE, and the simulation is equivalent to DNS. In the limit of infinite filter width, all the 
Reynolds stress is contained in sub-filter scales (i.e., in the residual stress), the resolved-scale 
Reynolds stress vanishes identically (because the field variables are time independent), the 
TFNSE reduce to the RANS equations (for stationary flow), and the simulation is equivalent 
to RANS, provided that the residual-stress model tends toward a Reynolds-stress model. In 
short, the larger A, the greater the burden placed on the model Mij ~ Rij, but also the greater 
the potential for grid coarsening. 

Finally, we consider the instantaneous turbulent kinetic energy, k(t), defined here as 


1 /‘ +1 

k(t) ■= - j <UiUi> XtV dz (22) 

The analogous computable quantity for TLES is the resolved turbulent kinetic energy, k(t ) 

l f +1 

k(t~j ~ J UjUi t>x,y dz 

Also of interest is the residual turbulent kinetic energy . kji(t) 

1 f +1 

= 2 J < ~ dba ' > x,y dz 


( 23 ) 


( 24 ) 
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For exact TLES, the instantaneous sum of resolved and residual turbulent kinetic energies 
is (only) approximately k(t). That is, 

k(t) k k(t ) + kn(t) (25) 

The time mean of the previous approximation, however, yields an exact equality; that is, 
k —< k(t ) >=< k(t ) > + < kn(t) >. 

V. DECONVOLUTION METHODS 

As suggested in the Introduction, deconvolution methods seek to approximate FUj pri- 
marily by mathematical means rather than to model by physical considerations. In 
this section, deconvolution methods are summarized briefly, the classical (spatial) ADM is 
presented, the ADM method is adapted for causal filters, and difficulties associated with 
causal filtering are addressed. 

A. Approximate Deconvolution 

A lucid review of deconvolution methods can be found in Adams and Stolz 2 , to which 
the reader is referred for details. Here, for completeness, we review the essentials. 

Let Uj-(£) be the temporal or (one-dimensional) spatial Fourier transform of the velocity 
field Uj, where £ = kA or £ = to A for spatial or temporal filters, respectively, and where k 
is the wavenumber and u is the circular frequency. Accordingly, A is either the spatial or 
temporal filter width. The Fourier transform of a filtered field is given by 

Ut) = H{Z)u j {0 (26) 

where H is the transfer function of the filter (e.g., Fig. 1). If H ^ 0, the filter is said to be 
invertible , in which case 

u j (0 = H- 1 (OUO ( 27 ) 

Using a standard trick from complex analysis, 

#~ 3 ~ ( ^_ j p j = 1 + (1- H) + (1 — H) 2 + ... + (1 — H) p + ... (28) 


14 


The geometric series on the right of Eq. 28 converges provided |1 — H\ < 1. Truncating the 
series at finite order p yields an approximate inverse of the transfer function, namely 

H~ l = 1 + (1 - H) + (1 - H ) 2 + ... + (1 - Hf (29) 

If H~ l were exact, then H' 1 * H = 1, and all scales would be faithfully recovered from the 
filtered fields. Exact inversion is not desirable for present purposes. Rather, one desires to 
recover the resolved scales 0 < |£j < £ c exactly, but to appropriately dissipate the energy 
in the unresolved scales £ c < j£| < n. Thus, for LES or TLES, the ideal is that H * H~ l 
approximates a sharp cutoff at £ c in Fourier space. 

Let F represent the filter operator in physical space, and let F k denote a /c- times iterated 
filtering operation. That is; 

Uj = Fuj 

Uj = FUj = F 2 Uj 


uf = F\ (30) 

By isometry between physical and Fourier space (Stolz and Adams 31 ), Eq. 29 can be ex- 
ploited to yield a deconvolution scheme in physical space, namely 


Uj « Vj = [/ + (/- F) + (/ - F) 2 + ... + (/ - FY)uj 

- t^r 

h=0 


(31) 


where Vj is the pth-order deconvolution approximation of Uj, and where the coefficients c* 
are reeudily determined to any order p by the binomial formula. Equation 31 represents 
but one of many possible approximate deconvolution techniques (See Adams and Stolz 2 ); 
however, it is advantageous in being linear and involving only multiply filtered quantities. 
For specificity, approximate deconvolutions of orders p = 0, p = 1, and p — 2 are given 
below. 


Vj = Uj (p = 0) 

v.j = 2 Uj - Uj (p = 1) (32) 

Vj = 3uj - 3 Uj + Uj (p = 2) 
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B. Spatial Approximate Deconvolution Method 


Stolz et al. 32 present two anisotropic spatial ADMs, which are reproduced below (in 
cdrrent parlance) for completeness. 


[Mi }ij = ViVj - UiUj 

(33) 

[M 2 }ij = vpTj - v^j 

(34) 


Model M 2 reduces to the scale similarity model (SSM) of Bardina 4 for p = 0 and accordingly 
is dubbed the generalized SSM (GSSM). Although Stolz and Adams find little difference 
between the two models in practical applications, only the latter is suited for TLES, for 
reasons addressed subsequently. 

It is well known that models of similarity type provide too little dissipation for practical 
applications to LES. The ADMs suffer from a similar deficiency without secondary regu- 
larization , , a type of high-order artificial dissipation that is accomplished by introducing a 
relaxation term to the right-hand sides of each momentum equation. In physical space, the 
relaxation term is of the form 

-X(uj ~ vj) (35) 

where x is an arbitrary non-negative parameter, to which the solution is not very sensitive 32 . 
In Fourier space, the relaxation term becomes 

— x(f - H * H~ l ) * Uj (36) 

Provided (7 — H * H ~ 1 ) is positive semi-definite, which is the case for symmetric spatial 
filters, the relaxation term is dissipative for all unresolved Fourier components. 

In the original spatial ADM, the cutoff is fixed at £ c = 2/3n. For this value, Stolz et al. 32 
report the ADM to give acceptable results for p = 3 and that deconvolution orders above 
5 do not significantly improve its performance. Consequently, they fix the deconvolution 
order at p = 5. Stolz et al. also observe that secondary regularization is equivalent to 
imposing a high-order filter upon the solution. According to Stolz et al. 32 , the effective 
order of the secondary filter is the product of the order of the primary filter and the order of 
the deconvolution plus one. For example, for the spatial ADM of Reference 32 , which exploits 
a 4th-order compact-difference filtering and p = 5, the secondary filter is formally of order 
24. Thus, in the classical ADM method, secondary regularization acts virtually as spectral 
filtering. 
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C. Temporal Approximate Deconvolution Method 


The proposed temporal ADM (TADM) is formally identical to the second of the ADMs 
proposed by Stolz et al.; that is, Eq. 34, where Vj is given by the truncated geometric series 
Eq. 31, except, of course, that now F denotes a causal filter operator. Because phase error 
introduced by temporal filtering results in the first and second terms of Eq. 33 being out of 
phase respective to one another, only the second ADM model (Eq. 34) is viable for present 
purposes (TLES). Henceforth, we refer symbolically to the residual stress modeled by the 
temporal version of [A^jy simply as My. The TADM also exploits secondary regularization 
formally identical to Eq. 35. 

Recall that Eq. 28, from which Eq. 31 is obtained, converges only for (1 — H\ <1. Fig 4 
shows this quantity for both the exponential filter of current interest and the Heaviside filter 
of Section II of Reference 27 . For the exponential filter |1 — H(£) | < 1 for all 0 < £ < oo. 
In contrast, this constraint is violated for high £ for the Heaviside filter. Consequently, the 
Heaviside filter is not a viable candidate for temporal deconvolution techniques of this type. 
Henceforth, we consider only the exponential filter (while recognizing that other filters yet 
to be examined could be superior in many respects). 

D. Adaptations of Original ADM for Temporal Deconvolution 

Straightforward adaptation of approximate deconvolution (Eq. 29) for the TADM is prob- 
lematic, chiefly due to the phase error generated by causal filtering. (These difficulties were 
partially masked in the Burger’s flow considered previously by the authors 27 , because of the 
intrinsic physical dissipation present at low Reynolds number.) For spatial filtering with a 
centered stencil, H, H~ l , and H * H~ l are each purely real, and the operator H * ET -3 , 
which is exploited for secondary regularization, is itself a low-pass filter. (See Fig. 1 of Stolz 
et al. 32 ) For causal filtering, each is complex; as a consequence, H * H~ l is poorly behaved. 
Figure 5 displays the moduli |H|, |H -3 |, and | H * H~ l \ obtained from the differential fil- 
ter (Eq. 8) and p-th order deconvolution (Eq. 31) for selected values of p. It is clear that 
the product | H * H~ l \ converges as expected toward the identity operator as p increases; 
however, convergence is slow, and for low-order deconvolution, conventional deconvolution 
actually increases energy in unresolved scales of motion. Moreover, as stated previously, for 
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deconvolution methods with applications to LES/TLES, one wants the product H * H~ 1 not 
to replicate the identity operator but to have the character of a high-order low-pass filter, 
so that large scales of motion are recovered faithfully even as small scales of motion are 
attenuated. 

These difficulties can largely be circumvented by defining the approximate deconvolution 
in the standard way as 

Vi = j2 c ^ +1) ( 3? ) 

k—0 

but with coefficients c* optimized to give both | H * H~ l | and \I — H * H~ l \ the following 
desirable properties: . 

1. *if- 1 (0)| = 1 

2. \H « 1 

3. | H*H~ l \ < 1 for all Q 

4. A minimal sum of absolute values of the low-order derivatives of | H * H~ l | at the 
origin 

5. \I-H*H~ l \ < 1 

The first constraint, imposed as strict equality, requires that the coefficients sum to unity; 
that is, Y7j=o c k — h The remaining constraints are “soft” in that they can be imposed in 
optimization packages via inequalities or minimization requests. The use of soft constraints 
prevents the system from being overconstrained whenever p is small. For example, for p = 3, 
a set of coefficients found by optimization strategies is 

{co, ci, c 2 , c 3 } = {0.99693, 0.503349, -0.684561, 0.184282} (38) 

Figure 6 shows the moduli \H\, |# _1 |, | H * H~ l \ and \I — H * H~ l \ for the 3rd-order 
deconvolution coefficients immediately above and r = 1. Similarly Fig. 7 displays the same 
information for r = 4. Note that in either case \H * H~ l \ has desirable low-pass properties. 
Specifically, it has a fiat plateau near the origin, so that large scales are recovered faithfully, 
yet small scales of motion are highly attenuated. 
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Figures. 6 and 7 also present the operator |1 — H * H~ } | which is germane to secondary 
regularization. That the moduli lie between zero and one implies that all scales (except 
the mean) are attenuated by the secondary regularization operator, with extreme damping 
for the higher scales. The ideal shape for this operator is an inverted low-pass filter, with 
a flat plateau near the origin. Whereas, the operator is indeed inverted low-pass, it has 
non-zero slope near Q. — 0, an undesirable trait. Thus, the TADM as presently configured 
may be more sensitive to the regularization parameter x than is the ADM method, which, as 
mentioned previously, is reported to be relatively insensitive 32 . In summary, the coefficients 
Cfc presented above should be considered sub-optimal. Optimized temporal deconvolution 
will (hopefully) be the subject of a follow-on paper. 


E. Governing System 


Recall that in the current time-filtered approach, the filter is imposed in differential form 
(Eq. 8). The governing equations thereby consist of evolution equations for the filtered 
velocities coupled to a set of evolution equations for additional filtered quantities as follows: 


dUj 

dij 

dui d(uiUj) 
dt dij 

Mij 


= 0 


dp 1 d 2 Ui d M i: 


dxi Redxjdij dij 


ViVj - ViVj 


= S^c fc u! 




—u {k+1) 
dt 1 

dvi 

~dt 


k = o 

- k -(fe+i) 

1 U ' (1 <k< p) 


A 

Vi - Vi 


(39) 


d , , ViVj - 


where the entire system is parameterized by the temporal filter- width A. As was shown in 
Pruett et al. 27 , the exact and modeled residual stresses vanish identically as A — * 0. Thus, 
as the filter width tends toward zero, the TFNSE tend smoothly toward the NSE, and TLES 
tends smoothly toward DNS. 
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The system above is presented without secondary regularization, which as mentioned 
previously, is a necessary form of artificial dissipation. The imposition of secondary regular- 
ization is simple; to the right-hand side of Eq. 39 expression 35 is appended. With secondary 
regularization, there are two arbitrary parameters, the deconvolution order p and the dis- 
sipation parameter x - The filter width A is not considered arbitrary; it can and should be 
established on the basis of physical considerations, as discussed later. 

At first inspection, it appears that the computational overhead for TLES is high, both 
in terms of storage and machine operations. In addition to the original storage required for 
the four primitive variables, storage locations are required for the nine (distinct) variables 
Vj and upv], as well as for the 3 p variables Uj k \ j — 2, ...,p+ 1. Thus, for p = 0 and p = 3, for 
example, the total storage for TLES is about three times and seven times, respectively, that 
of a simulation without the model. However, if TLES permits significant grid coarsening 
relative to fully resolved DNS, then the net storage savings could be significant. For example, 
for p = 3, a TLES coarsened by a factor of four in each spatial dimension and time relative 
to DNS would necessitate only 10 percent of the memory required by the DNS. TLES is 
more benign with regard to computational effort. Although the governing system has grown 
from three evolution equations to 9 + Zp equations (vj do not involve evolution equations), 
the additional equations are linear and involve no spatial differentiation. Hence, even for 
p = 3, the computational cost of evaluating all additional filter equations remains less than 
that of evaluating the three momentum equations. Further reductions in computational 
overhead, for example, would be possible if one were willing to settle for an isotropic residual- 
stress model. No attempt to optimize either storage or computational effort has been yet 
undertaken, and it is likely that the TADM can be streamlined considerably. 

VI. REFERENCE DNS 

For several reasons, plane-channel flow affords an ideal test case against which to bench- 
mark present methodology. Channel flow' involves wall effects, which are highly anisotropic, 
yet lends itself to fully spectral numerical methods, highly prized because of their accuracy 
and computational efficiency. Moreover, in its stationary turbulent state, channel flow is 
homogeneous in two of three spatial dimensions, which reduces the data output of post- 
processing. For these reasons, channel flow has been studied extensively 14-17,34 , and a highly 
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accurate data base is available (Moser et. al 17 ). By convention, we will refer to this data 
base as KMM, as it originated with the work of Kim et al. 15 As a point of reference, Fig. 8 
presents the streamwise velocity (u } ) profiles for the KMM cases of Re T = 180 and Re T = 590 
vs. that of the laminar state, which is parabolic in the transverse coordinate. Here, Re r is 
the Reynolds number based upon the mean friction velocity u T (to be defined shortly), and 
channel half width. To ensure consistency of results, we have repeated the DNS of Moser et 
al. for nominal Re r = 180. 

A. Numerical Method 

The DNS is conducted by solving the NSE (that is, no residual-stress model) with the 
spectrally accurate channel-flow code TRANSIT of Gilbert and Kleiser 12 . The Navier-Stokes 
equations are solved in spectral space using Fourier expansions in the periodic streamwise 
(x) and transverse ( y ) directions and Chebyshev polynomial expansions in the wall-normal 
( z ) direction. Collocation points are equally spaced in the periodic dimensions and are 
distributed at the Gauss- Lobatto points in the wall-normal direction. Nonlinear products 
are computed pseudospectrally by fast transform methods. An option exists in the code to 
run with dealiasing according to the 3/2 rule or without explicit dealiasing. Another option 
exists regarding the assumption of transverse symmetry; per-step computational efficiency 
is improved by a factor of two with symmetry enforced. Time advancement is semi- implicit 
and couples Crank-Nicolson splitting of diffusion terms with third-order low-storage Runge- 
Kutta time advancement of advection terms. A hallmark and innovation of TRANSIT 
is the enforcement of the divergence-free condition to machine (double) precision through 
an influence-matrix technique. The computational domain is a box of dimensions Li = 
2 n / ai , Z /2 = 27t/q; 2) and L 3 = 2 in the streamwise, transverse, and wall-normal directions, 
respectively. TRANSIT, which performs temporal DNS, permits the computational box to 
convect at phase velocity CPH. All computations discussed herein, both DNS and TLES, 
were performed with CPH=0.0, for reasons explained subsequently. 

For some simulations, the flow is initialized to a perturbed laminar state, for which the 
streamwise velocity profile is parabolic (Fig. 8) with centerline value unity; that is, u(z) = 

1 — z 2 . For others, to spare computational effort, the simulation is initialized to a nominally 
turbulent state. In either case, the flow is allowed to evolve through an initially transient 
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period to stationarity, at which point statistical post-processing begins. Which of these 
initial conditions is exploited for a particular simulation will be reported in context. Without 
the enforcement of transverse symmetry, randomness in the initial conditions results in 
slightly asymmetric velocity and Reynolds- stress profiles. Rather than to enforce symmetry 
during the simulation, we prefer to allow asymmetry but to compute statistical quantities by 
averaging the upper and lower channel halves, which can be considered different realizations. 
Thus, the additional computation time per step is offset by the need for fewer time steps. 

Channel-flow simulations are non-dimensionalized in several ways. To avoid confusion, 
we draw attention to three different methods for expressing the Reynolds number. In each 
case, Re = U6/ u, where U is the reference velocity, 5 is the channel half width (assumed 
■without loss of generality to be unity), and v is the kinematic viscosity. The specifications 
differ in their reference velocities U. itej am , Re b) pp, and Re T result from using the laminar 
centerline velocity C/| am , the laminar bulk velocity N bubci and the friction velocity u T , 
respectively, as the reference velocity, where u T = \J ^ d< ^ > ■ Without loss of generality, we 
presume I7 lam = 1, in which case {7 bulic = 2/3; hence, Re bulk = 2/3Re lam . The nominal 
Re r = 180 case of Moser et al. corresponds to Re bu j b = 2800 (Rejg^ = 4200). The flow 
is driven either by a stream-wise pressure gradient or by requiring the mass-flow rate to 
be constant. The latter option is preferred computationally, because the flow settles more 
quickly to stationarity (N. Adams, personal communication). Moreover, for this option, the 
mass-flow rate remains 4/3 at all times, as can be readily derived by integrating the laminar 
initial profile between -1 and +1 with I/} am = 1. K no-wing the mass-flow rate in advance is 
advantageous in ensuring consistency among competing scalings. 

Finally, we note that Re r is the time mean of the instantaneous quantity Re T = u T 5/v , 
where u T is the instantaneous friction velocity, defined asu 7 = y v d<U g> x ' y ■ 

B. DNS Test Cases 

Parameter values for the DNS test cases are summarized below in Table I. Here N x , 
N v , and N z are the numbers of collocation points in the respective coordinate dimensions. 
Note that for TRANSIT, Rej am = 3/2Re bu j b is a fixed input parameter, but that Re T is a. 
time-dependent output result of the simulation, whose long-time mean is Re T . 
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TABLE I: Parameter values for reference DNS cases. 


Case 

N x x N y x N z 

At 

C*j 

a 2 

Re buUc 

Re r 

dealiased? 

KMM 

128 x 128 x 129 

— 

0.5 

1.5 

2800 

178.1 

yes 

DNSa 

128 x 128 x 129 

0.01 

0.5 

1.5 

2800 

176.2 

yes 

DNSb 

128 x 128 x 129 

0.01 

0.5 

1.5 

2800 

« 162 

no 


C. Results 

To provide a well-documented reference flow against which to compare TLES results, we 
discuss the DNS results, with close attention to stationarity, homogeneity, and anisotropy, in 
addition to presentations of mean velocity profiles and mean Reynolds stresses. For brevity, 
we refer to simulations with resolution 128 x 128 x 129 as having resolution 128 3 . 

1. Homogeneity 

In temporal DNS of plane channel flow, the flow is (by definition) periodic and homoge- 
neous in the streamwise dimension. Channel flow is also periodic in the spanwise dimension; 
however, homogeneity must be established. Figure 9 displays < ui > t .x (y) at two (ar- 
bitrary) wall-normal stations: z = 0 and 2 = .9988. The spanwise variations are shown 
relative to their spanwise mean values. Temporal window averages are obtained for A = 50 
and A = 200. As the window size increases the amplitude of the spanwise deviation about 
the mean appears to diminish, which suggests that the flow is homogeneous in y for an ap- 
propriately long time scale (somewhat longer than present computational limitations allow) . 
We now exploit homogeneity to reduce the turbulent statistics to functions of one dimension 
(2) by averaging spatially over wall-parallel planes. 

2. Stationarity 

Figure 10 displays the instantaneous and mean values of two turbulence statistics of 
interest: Re r and turbulent kinetic energy k for the time interval 0 < t < 200. These results 
suggest that the flow is stationary; however, rigorous confirmation of stationarity would 
require a time interval considerably longer than what practical computational limitations 
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allow. 


S. Aliasing Error Control 

An initial surprise was the effect of aliasing error upon the flow statistics. Cases DNSa and 
DNSb differ only in the parameter that invokes or bypasses aliasing error control. However, 
as shown in Fig. 11, their results differ by approximately 16 in Re T . Moreover, such a 
difference in Re T induces noticeable differences in mean velocity profiles. The significant 
differences between the aliased and dealiased results suggest that spectral simulations at 
Re r = 180 of 128 s resolution are adequately resolved but not highly resolved. Given that 
computational limitations typically force simulations to be conducted at marginal resolution, 
it is imperative that aliasing error be controlled explicitly. 

For TRANSIT, dealiasing is performed unconventionally by increasing the resolution by 
a factor of 3/2 in each dimension prior to the transform to physical space, where nonlinear 
products are computed. Following the evaluation of nonlinear terms, resolution is returned 
to its nominal values by spectral truncation in Fourier space. This approach, which saves 
storage, is costly in another respect, namely that the per-step computational time nearly 
quadruples with explicit dealiasing. The necessity of controlling aliasing error also has 
connotations for TLES, which are addressed later. 

4- Mean Profiles 

Figure 12 compares the mean streamwise velocity and its wall- normal derivative for the 
present simulation against that of the standard of Moser et al. The streamwise velocity 
profile is presented both in standard units and in wall units, uf [= u-l/ut ] vs. z + [= 
zu T /u\. Whereas the presentation in standard units requires rescaling the KMM data, the 
presentation in wall units requires rescaling current results. That the results are virtually 
indistinguishable ensures not only that the two independently developed codes are giving 
consistent results, but also that the different scaling conventions have been interpreted in a 
mutually consistent manner. 
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5. Reynolds Stresses and Anisotropy 


Figure 13 compares mean Reynolds stresses from the present simulation with the standard 
of Moser et al. for nominal Re T = 180. The agreement is very good, which validates both 
the algorithm and the parameter settings. 

The trace elements of the mean Reynolds-stress tensor afford a measure of anisotropy, 
which we take to be the ratio of twice the streamwise component to the sum of the remaining 
diagonal components, the value of which is presented in Fig. 14 as a function of wall-normal 
coordinate. For an isotropic flow, this ratio is identically one. As the figure shows, channel 
flow is highly anisotropic near the walls, being strongly dominated there by streamwise 
velocity fluctuations. 

In Eq. (7.10) of Pope 22 , the author presents a shear stress balance for turbulent channel 
flow. In dimensionless form with current scalings and (without loss of generality) p — 1, the 
total shear stress r depends only on 2 as follows: 


T(z) _ _ 713 

u 2 dz u 2 


(40) 


where r 13 =< u[u' 3 > is the (1,3) component of the Reynolds-stress tensor. Moreover, in 
Eq. (7.13), Pope shows the total stress profile to be linear. That is, the sum of viscous and 
Reynolds stresses across the channel is linear. Figure 15 confirms this relationship for case 
DNSa. 


VII. RESULTS OF TLES 

The current section consists of three subsections. In the first, three coarse-grid simulations 
of channel flow at nominal Re T = 180 are presented to establish a baseline for evaluating 
the effectiveness of TLES with the TADM, the results of which are presented in the second 
subsection for Re T = 180. The final subsection presents TLES results of channel flow for 
nominal Re r = 590. For each case, results are referenced to those of DNS at equivalent 
Re r , as discussed in the previous section. All simulations were conducted with CPH=0.0, 
a necessity, because non- zero phase velocity produces a Doppler shift in the effective filter 
frequency 23 that interferes with the interpretation of some results. 
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A. Baseline Cases 


Parameter values for the baseline coarse-grid cases at nominal Re T = 180 are presented 
in Table II below. For all cases a?i = 0.5, a 2 = 1.5, and = 2800 (-^ e iaxn = 4200), in 

accordance with the values for the KMM case of (nominal) Re r = 180. Whenever temporal 
filtering is involvee, the filter- width ratio r = 8. The reader is reminded that, whereas 
(P^lam) is an input parameter, Re r and its mean Rtr are output values. The baseline cases 
are designed not to assess the effectiveness of the TADM, but rather to establish coarse-grid 
reference simulations against which the method may be evaluated. For all baseline cases, 
the initial condition is derived from fine-grid DNS results, which are interpolated onto the 
coarser grid. As such, the initial condition induces transients and does not identically satisfy 
zero mass divergence. The latter is satisfied to machine (double) precision at the end of the 
initial time step; however, it takes quite some time for transients to settle. The baseline 
cases represent simulations with no model (BASEa) and simulations with a minimal model 
(BASEbc), namely, a temporal scale- similarity model (TSSM), which is the degenerate case 
of the TADM for p = 0 27 . The TSSM is the temporal analog of the Bardina model 4 , which 
is known to be insufficiently dissipative. Cases BASEbc differ only in the de-aliasing option. 

In Stolz et al. 32 , the authors reported that, prior to the advent of the ADM, of all 
residual-stress models tested for channel flow, the most effective was no model whatsoever. 
This situation was improved substantially by the introduction of the ADM, which in our 
judgment, now represents the de facto state-of-the-art in residual-stress modeling, certainly, 
at least, for channel flow. Consequently, Table II also includes parameter data for Case 
ADM-180a of Stolz et al. 32 , which we refer to as SAKl80a. The baseline simulations and that 
of SAK represent a suite of coarse-grid simulations against which to evaluate TLES/TADM. 
More specifically, TLES/TADM should be considered a moderately effective methodology 
if it improves substantially upon the “no model” case (BASEa), and very effective if it 
produces results roughly equivalent in accuracy to the standard (spatial) ADM method at 
commensurate resolution. 
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TABLE II: Parameter values for baseline coarse-grid simulatons. 


Case 

N x x N y x N z 

At 

# e bulk 

Re T 

dealiased? 

model 

(dp) 

BASEa 

32 x 32 x 33 

0.04 

2800 

203 

yes 

none 

NA 

BASEb 

32 x 32 x 33 

0.04 

2800 

182 

no 

TSSM 

(8,0) 

BASEc 

32 x 32 x 33 

0.04 

2800 

203 

yes 

TSSM 

(8,0) 

SAK180a 

32 x 32 x 33 

— 

2800 

173 

yes 

ADM 

NA 


1. Case BASEa: No Residual-Stress Model ; with Dealiasing 

Figure 16 displays computed Re r and turbulent kinetic energy k for Case BASEa. This is 
also referred to as the “no-model” TLES. which is essentially coarse-grid DNS with dealias- 
ing. The figure suggests that the flow is stationary by t — 100: consequently mean quantities 
presented below are computed by temporal averages taken over the interval 100 < t < 500, 
as shown by the length of the horizontal line in Fig. 16 and subsequent figures. As discussed 
previously, mean data are computed by also averaging over homogeneous dimensions and 
exploiting symmetry. 

Figure 17 presents the mean streamwise velocity and its wall-normal derivative vs. wall- 
normal coordinate z for Case BASEa. Although the “no model” simulation appears to 
produce a nearly correct mean velocity according to Fig. 17(a), the slope of the velocity- 
profile at the wall is approximately 30 percent in error, as indicated in Fig. 17(b)-(c). The 
error in wall velocity gradient results in Re r = 203.2 for Case BASEa vs. Re r = 178.1 
for the DNS. Figure 17(c) presents the same information as Fig. 17(a), except expressed 
in wall units. More precisely. In Fig. 17(c), each profile is scaled according to the relevant 
calculated friction velocity u r = i?e T //?e] am : that is, according to Re T = 178.1 for DNS and 
Re T = 203.2 for BASEa. This convention has the effect of shifting the origin so that profiles 
match at the wall but deviate due to accumulated error away from the wall. 

Figure 18 presents the four principal components of the Reynolds-stress tensor for Case 
BASEa, relative to the Reynolds stresses of the reference DNS of KMM. In Stolz et a!. 32 , 
a precedent was established for displaying trace elements of the Reynolds-stress tensor as 
fluctuation velocities R\( 2 ■ We adopt the same convention in Fig. 18 and subsequently. 
Relative to DNS, the “no model” case over-predicts all components of Ty, with the (2,2) 


27 



component some 60 percent in error. Here, to compare Reynolds stresses in absolute terms, 
both the DNS and BASEa values have been scaled by the same reference velocity, namely 
u T of the DNS. The customary rescaling of Reynolds stresses by the actual u T (in this Case 
xir = Re T /i?ej am = 203.2/4200) gives the appearance of a better fit of the model to the 
DNS reference data, but, in our view, is misleading because it masks the effects of under- 
or over-estimating the wall stress t w , 

2. Case BASEb: TSSM Residual-Stress Model without Dealiasing 

Figures 19, 20, and 21 display the same information as in the previous case, but for 
Case BASEb. Figure 19 suggests that the flow is somewhat transient until t « 200, but is 
essentially stationary for t > 200. Average Re r for 200 < t < 500 is 181.9. To approximate 
k(t) for Case BASEb, we have exploited Eq. 25, albeit with modeled residual stress Mu in 
lieu of Ra- 

For Fig. 21, Reynolds stresses for Case BASEb are approximated as the sum of the 
resolved-scale Reynolds stress and the mean modeled residual stress, according to Eq. 21 
with Mij replacing Rj. The contribution of the modeled residual stress to the total Reynolds 
stress is also presented. Reynolds stresses are considerably in error near the walls, where T n 
is under-predicted and r 2 2 and r 33 are over- predicted. In particular, the (2,2) component is 
over-predicted by about 60 percent. Errors in Reynolds stress induce considerable mean-flow 
error, as shown in Fig. 20. 

From these results and the results of the DNS of the previous section, we conclude that 
dealiasing is essential. 

3. Case BASEc: TSSM Residual-Stress Model; with Dealiasing 

Figures 22, 23, and 24 display the same information as previously, but for Case BASEc, 
for which the model is the TSSM. In general, it appears that the p = 0 case performs no 
better and no worse than the no-model case. That all three coarse-grid baseline test cases 
predict Reynolds stresses and turbulent kinetic energy to be in excess of that of the DNS 
(Fig. 10) confirms that the baseline models are insufficiently dissipative as expected. 
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B. TLES with Pull TADM for Nominal Re T = 180 


Here we present results of TLES with a full TADM for nominal Re T = 180; that is, 
the residual-stress “model’’' incorporates both high-order deconvolution and secondary reg- 
ularization. For each case the deconvolution is of order p = 3 with “tuned” coefficients, 
as discussed in Section V. The parameters of the test cases are given in Table III below. 
As for the results of the previous section, ai = 0.5 and 0:2 = 1.5 for all test cases, and 
■^ e bulk = 2800 (-ftejam = 4200), for congruence with the results of KMM at Re r = 180. 
The filter-width ratio r and regularization coefficient x are provided in the last column of 
the table. Unlike the baseline cases established above, the initial condition for all TLES 
cases at Re T = 180 is that of a randomly perturbed laminar channel flow. 


TABLE III: Parameter values for TLES cases. 


Case 

N x x N y x N z 

At 

^bulk 

Rz t 

(r,x) 

TLES 180a 

32 x 32 x 33 

0.04 

2800 

172.5 

(08,0.5) 

TLES 180b 

32 x 32 x 33 

0.04 

2800 

179.2 

(08,0.3) 

TLES 180c 

32 x 32 x 33 

0.04 

2800 

184.7 

(08,0.2) 

TLESl80d 

32 x 32 x 33 

0.04 

2800 

162.0 

(16,0.1) 


For brevity we present graphically the results of Cases TLESb, TLESc, and TLESd only 
(Figs. 25-34). All TLES results except TLESd show marked improvement relative to the 
corresponding results of the coarse-grid baseline Cases BASEa-c. In particular, the TLES 
meanflow profiles very nearly match those of the DNS. Secondary regularization provides 
an adequate sink for energy at marginally resolved scales, as evidenced by decreases in 
turbulent kinetic energy k and Re T relative to the baseline cases. At present, the value of 
the regularization parameter x is arbitrary, and an appropriate value is found by numerical 
experimentation (as is also true for the conventional ADM). For fixed r, variations of x by a 
factor of 2.5 do not change the results dramatically, as can be seen by comparing the x = 0.5 
and x = 0-2 cases for r = 8. The insensitivity to x of the conventional (spatial) ADM is 
discussed in Stolz et al. 32 . As expected, the current TADM appears to be more sensitive to 
X than the ADM. Acceptable results, however, are obtained for a range of values of y. For 
Case TLESd, the level of dissipation imposed by the large filter width and the secondary 
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regularization is quite high, and it is clear that the combination is overly dissipative. 

In contrast to the ADM of 32 , for which the filter width is fixed, the current TADM 
is parameterized by the filter-width ratio r. As implied in Figs. 6 and 7, the secondary 
regularization operator becomes more dissipative with increasing r at roughly a linear rate. 
Improvement of the secondary regularization operator for the TADM will be the subject of 
follow-on research. 

Figure 28 below establishes for Case TLESc that TLES of channel flow follows a stress 
balance analogous to that shown in Fig. 15 for Case DNSa, except that there is an additional 
contribution to the total mean shear stress through the mean modeled residual stress. Here, 
the velocity scale is the computed (rather than nominal) value of u T = 179.2/4200. Low-level 
oscillations in the linear profile suggest, however, that 32 s spatial resolution is marginal. 

In contrast to the baseline cases, all TLES cases for Re r = 180 axe initialized to a lami- 
nar state subject to random perturbations. Thus, the computation proceeds from perturbed 
laminar flow, through transition, to a stationary turbulent state. Transitional flow is char- 
acterized by a strong spike in turbulent kinetic energy k , as evidenced in Figs. 25, 29, and 
32. In the initial stages of transition, the flow is primarily two-dimensional. Late-stage tran- 
sition is characterized by loss of two-dimensionality and the prevalence of three-dimensional 
structures. In the final breakdown stage, k “relaxes” toward a lower value characteristic of 
a fully turbulent state. 

In general, it is difficult for LES to properly treat transitional flow. On the one hand, 
the standard Smagorinsky model is excessively dissipative and typically inhibits transition 
altogther. On the other hand less dissipative methods tend to blow up during the spike in k 
associated with transition. That TLES, properly tuned, can allow for transition, yet settle 
into a statistically nearly-correct turbulent state, is noteworthy. 

C. TLES with Full TADM for Nominal Re T = 590 

Channel flow at nominal Re T = 180 is generally considered “barely turbulent.” Conse- 
quently, it is highly desirable to validate TLES at much higher Reynolds number, for which 
we consider nominal Re r = 590. Table IV summarizes the parameter values of the present 
TLES relative to those of the reference DNS of KMM for nominal Re r = 590. The table 
also presents parameter data regarding channel-flow DNS and LES conducted at the same 
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nominal Reynolds number by Stolz, Adams, and Kleiser (SAK) 32 . More precisely, all sim- 
ulations were conducted at = 10935 (Rej am = 16400) with box parameters au = 1 

and 02 = 2. For this parameter set, KMM obtained Re T = 587. Independent replication of 
the DNS by SAK, who used TRANSIT and the same spatial resolution as KMM, resulted 
in Re T = 586. Thus, very high confidence exists in the reference DNS data. 

Course-grid simulation with no model (other than explicit dealiasing) by SAK yielded 
Re r = 633, an error by about eight percent, which corresponds to an overshoot of approxi- 
mately 16 percent in the mean wall-normal derivative of the streamwise velocity at the wall. 
In constrast, SAK obtained Re T = 574 and Re T = 587 for their LES cases 590a and 590b, 
errors of just over two percent and zero percent, respectively. Cases SAK590a and SAK590b 
each exploit the standard (spatial) ADM and differ only in resolution. 


TABLE IV: Parameter values for TLES and reference simulations at nominal Re r = 590. 


Case 

N x x Ny X N z 

Re T 

model 

DNS (KMM) 

384 x 384 x 257 

587 

NA 

DNS (SAK) 

384 x 384 x 257 

586 

NA 

LES (SAK) 

48 x 64 x 65 

633 

none 

LES (SAK590a) 

48 x 64 x 65 

574 

ADM 

LES (SAK590b) 

72 x 96 x 97 

587 

ADM 

TLES590 

48 x 64 x 65 

564 

TADM 


We consider TLES at the same resolution as SAK590a with the tuned (p — 3) TADM 
coefficients of Eq. 38 and parameter values At = 0.04, r = 8, and x = 0.3. Figure 35 presents 
the evolution of instantaneous Re T and k for Case TLES590, from a nominally turbulent 
initial state. Figures 36 and 37 present the mean streamwise velocity (and its wall-normal 
derivative) and the Reynolds- stress distributions, respectively, for Case TLES590. It is 
fair to say that TLES with the TADM performs nearly as well as the standard ADM at 
commensurate resolution. It is worth noting that, in contrast to the TLES cases at nominal 
Re r = 180, no numerical experimentation was performed to optimize parameters. That 
is, optimal values obtained for Re r = 180 were simply translated directly to R.e r = 590. 
However, at this Reynolds number and with these parameter values, the simulation had 
insufficient dissipation to survive transition. Mean quantities are computed by averaging 
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over 500 < t < 2000, following the settling of initial transients. 


VIII. DISCUSSION 


The computational sayings of TLES relative to DNS are significant. For example, the 32 3 
TLES cases at nominal Re T = 180 ran to completion in approximately one one-hundredth of 
the time required for DNS at a resolution of 128 3 . The speedup comes not only because of 
grid coarsening by a factor of four in each spatial dimension, but also by similar coarsening 
in time, for an aggregate workload reduction factor of 4 4 = 256. However, at commensurate 
resolutions, the computational overhead of TLES relative to DNS is approximately a factor 
of two, for a net gain in efficiency of roughly 10 2 . At Re r = 590, the computational gains 
of TLES are even more significant, as can be seen by comparing resolutions for DNS and 
TLES in Table IV. 

It is useful to establish a practical rule of thumb in regard to the choice of r (which 
effectively reduces the number of model parameters). One approach is to consider the 
ratio of the filter width A to the temporal integral scale f, a scale characteristic of the 
turnover time of the largest eddies. Figure 38 presents selected mid-plane (z = 0) time 
traces of the velocity components from Case TLES590. Associated with any time series is 
its autocorrelation function p(s) (Fig. 39), which in general, is given by by 


P(s) 


< u(t)u(t + s) > 
< u(t) 2 > 


(41) 


Several measures of the integral scale f can be extracted from the autocorrelation function. 
The favored is 


r oo 

f = / p(s)ds (42) 

Jo 

provided that the integral converges. If the integral does not converge, it is sometimes 
useful to let the first zero crossing of the autocorrelation function represent f. Figure 39 
presents a composite autocorrelation function derived by averaging the N x * N y individual 
autocorrelations, each associated with a time series of Uj recorded at one of the mid-plane 
grid points. A time-series is displayed for ui in Fig. 38. The autocorrelation functions of the 
velocities are quasi-periodic in time, as is to be expected. Although strict temporal period- 
icity is destroyed by the nonlinearity of the governing system, enforced spatial periodicity 
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nonetheless induces temporal quasi-periodicity. Because of the large mean streamwise com- 
ponent, Eq. 42 fails to converge for uj. Consequently, we use the first zero crossings from 
the autocorrelations of U 2 and U 3 as representative of the integral scale; thus, t % 2. The 
temporal resolution of Case TLES590 is such that there are approximately 50 time steps 
per laxge-eddy turnover. Consequently, r = 8 translates to a filter width of approximately 
16 percent of the integral scale. 

This suggests a possible overarching rule of thumb. TLES should behave as RANS when 
the filter width A is many times the integral scale (provided the SFS model is adequate at 
large A). For A a sizable fraction of f , say 50 percent, then TLES should effectively act as 
unsteady RANS (URANS). At the other extreme, for A quite small, say 1-2 percent of f, 
TLES is virtually equivalent to DNS. Finally, for A values that range from, say, 5-25 percent 
of f , TLES is the temporal analog of LES. 

Equivalently, the distinctions among DNS, TLES, URANS, and RANS, can be drawn 
based upon the relative contributions of the residual stress to the Reynolds stress. As dis- 
cussed in Section IV for RANS, the residual stress contributes 100 percent to the Reynolds 
stress. For DNS, the residual stress vanishes, and thus contributes nothing to the Reynolds 
stress. Between these extremes lie TLES and URANS. We would suggest that contributions 
of 5-50 percent of the residual stress to the Reynolds stress define TLES, and contribu- 
tions greater than, say, 50 percent define URANS. These domain boundaries are admittedly 
arbitrary and are offered simply as talking points. 

Finally, a few comments regarding the present SFS model are in order. There are two com- 
ponents to the model: the residual-stress model and secondary regularization. Experience to 
date suggests that the temporal residual-stress model is more than adequate. Deconvolution 
orders of 2, 3, and 5 have been examined, as have a variety of deconvolution coefficients and 
filter widths. In all cases the instantaneous and mean residual stresses are well-defined, with 
mean components that have qualitatively correct profiles. In sum, temporal deconvolution 
appears to admit a family of robust residual-stress models. 

In contrast, temporal secondary regularization remains problematic, apparently due to 
the phase lag associated with causal filtering. As artificial damping (secondary regulariza- 
tion) is an essential component of the present model (or of any generalized similarity model), 
the issue requires further study and is the subject of an addendum to this report entitled 
“Optimized Temporal Deconvolution.” 
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IX. CONCLUSIONS AND FUTURE WORK 


The present paper establishes temporal LES (TLES) as a viable alternative to conven- 
tional (spatial) LES. The methodology exploits a causal time-domain filter expressed in 
differential form, which yields a governing system of equations explicitly parameterized in 
terms of the temporal filter width A. Sub- filter scales (SFS) are approximated by a temporal 
variant (TADM) of the approximate deconvolution model (ADM) of Stoltz and Adams 31 . 
The method is demonstrated by simulating plane-channel flows at Re r = 180 and Re r = 590. 
Results axe compared against the benchmark channel-flow results of Kim et al. 17 , obtained 
by well-resolved DNS with spectral methods. 

Although there is likely considerable room for improvement in the current TADM, present 
results of TLES /TADM are nearly as good for channel flow as those obtained by LES/ ADM 
and reported in Stolz et al. 32 (for commensurate resolution). 

TLES enjoys some conceptual and practical advantages, which could make it the method 
of choice in certain applications. First, through the parameter A, TLES provides a formal 
link among DNS, URANS, and RANS methodologies. Second, whereas spatial filtering is 
problematic for unstructured grids because of wide variations in grid size, temporal filtering, 
particularly in differential form, is independent of spatial grid resolution. Finally, present 
methodology provides direct means for comparing results of TLES with DNS or experiment, 
obviating the need for a posteriori analyses. 
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FIG. 3: Transfer functions of continuous and high-order discrete exponential filters. 
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FIG. 4: Quantity |1 - H\, which determines convergence of Eq. 28 for exponential and Heaviside 
filters showing that Heaviside filter is not suitable for deconvolution. 
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FIG. 7: Moduli of transfer function of exponential filter (|fT|), its approximate inverse (| H -1 1), their 
product (\H and secondary regularization operator (|1 — H *H~ l \) for tuned deconvolution 
of order p = 3 and filter-width ratio r = 4. 
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FIG. 8: Mean streamwise velocity profiles of KMM vs. parabolic profile of laminar state. (KMM 
data are replicated symmetrically about centerline for clarity.) 
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FIG- 9: Spanwise variation in average streamwise velocity < uj > £jX (y), for temporal windows of 
A = 50 and A = 200. Amplitude of spanwise variation diminishes with increasing window size, 
suggesting spanwise homogeneity on a sufficiently large time scale. LEGEND: solid line (A = 200); 
dashed line (A = 50); dashed-dotted line (mean). Upper figure: z — 0.0. Lower figure: 2 = 0.9988. 
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FIG. 12: Mean streamwise velocity < u\ > and its wall-normal derivative in standard and wall 
units for Case DNSa relative to KMM reference data. LEGEND: solid line (KMM); dashed line 
(DNSa). (KMM data axe replicated symmetrically about centerline for clarity.) 
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nal Rer ~ 180. LEGEND: solid line (KMM); dashed line (DNSa). (KMM Data are replicated 
symmetrically about centerline for clarity.) 
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FIG/ 15 : Shear stress balance for Case DNSa. following Pope 22 . Figure confirms theoretical result 
that sum of viscous and Reynolds stresses is linear for turbulent channel flow. LEGEND: solid line 
(total shear stress [= r/u 2 ]); dashed line (viscous stress [= duj / dz ]); dashed-dotted line (Reynolds 
stress [=Ti3/u 2 ]). 
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FIG- 17: Mean streamwise velocity profile < u\ > vs. z and its wall-normal derivative for “no 
model” coarse-grid Case BASEa relative to reference data of Moser et aL, presented in both stan- 
dard and wall units. LEGEND: solid line (KMM); dashed line (BASEa); symbols (BASEa, wall 
units) . 
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FIG. 18: Fluctuation velocities and Reynolds stresses vs. z for Case BASEa relative to reference 
DNS data of Moser et al.: solid line (reference DNS), dotted line (BASEa scaled absolutely by 
Ur — Rer/Re lam = 178/4200 of reference DNS); dashed line (BASEa scaled relatively by derived 
Ur — Re T /Rei Qin ~ 203/4200). 
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FIG. 20: Mean streamwise velocity profile < u\ > vs. z and its wall- normal derivative for Case 
BASEb relative to reference DNS data of Moser et al. LEGEND: solid line (KMM); dashed line 
(BASEb); symbols (BASEb, wall units). 





















FIG- 23: Mean streamwi se velocity profile < u\ > vs. z and its wall-normal derivative for Case 
BASEc relative to reference DNS data of Moser et al. LEGEND: solid line (KMM); dashed line 
(BASEe); symbols (BASEc, wall units). 
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FIG. 26: Mean streamwise velocity profile < ui > vs. z and its wall-normal derivative for Case 
TLESl80b relative to reference DNS data of Moser et al. LEGEND: solid line (KMM) ; dashed line 
(TLES180b); symbols (TLESl80b, wall units). 
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FIG. 28: Shear stress balance for Case TLESl80b, following Pope 22 . Figure confirms theoretical 
result that sum of viscous and Reynolds stresses is linear for turbulent channel flow. LEGEND: 
solid line (total shear stress [= r/u 2 )]); dashed line (viscous stress [— 1/Re r du^ /dz))\ dashed- 
dotted line (resolved Reynolds stress [= fi 3 /t£ 2 ]); dotted line (mean modeled residual stress [= 

- < M 13 > /u?j). 
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FIG. 29: Evolution of Re r and k for Case TLESl80c. LEGEND: solid line (instantaneous); dashed 
line (mean); dashed-dotted line (fc); dotted line (kji). 




FIG. 30: Mean streamwise velocity profile < u\ > vs. z and its wall- normal derivative for Case 
TLESl80c relative to reference DNS data of Moser et al. LEGEND: solid line (KMM); dashed line 
(TLES180c); symbols (TLESl80c. wall units). 
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FIG. 31: Fluctuation velocities and Reynolds stresses vs. z for Case TLESl80c relative to reference 
DNS data of Moser et al. LEGEND: solid line (reference DNS); dotted line (< Afy > -hfy , TLES); 
dashed line (contribution of < Mij >, TLES). 


FIG. 32: Evolution of Re r and k for Case TLESd. LEGEND: solid line (instantaneous); dashed 
line (mean); dashed-dotted line ( k)\ dotted line (kn). 











FIG. 33: Mean stream-wise velocity profile < m > vs. z and its wall-normal derivative for Case 
TLESl80d relative to reference DNS data of Moser et al. Symbols denote values at collocation 
points. 
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FIG. 36: Mean streamwise velocity profile < uj > vs. z and its wall-normal derivative for Case 
TLES590 relative to reference DNS data of Moser et aL LEGEND: solid line (KMM); dashed line 
(TLES590); symbols (TLES590, wall units). 
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FIG- 38: Mid-plane (z — 0) time-trace segments of velocity components for Case TLES590. LEG- 
END: upper solid line (u x ); fourth (u[ 2) - 0.15); third (uf } - 0.3); second (u 2 4* 0.15); lower 
(U 3 ). Selected curves shifted for clarity. Upper three traces collectively show smoothing effects of 
successively filtering with causal exponential filter. 
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1 Introduction 

Recently, Stolz and Adams (1999) revealed an approximate deconvolution model (ADM) (see 
also Adams and Stolz (2001)) for large-eddy simulation (LES), which has performed well in 
simulating flows as diverse as incompressible channel flow (Stolz et al (1999a)) and supersonic 
compression-ramp flow (Stolz et al (1999b); Adams (2002)), and appears to represent the 
state-of-the-art in residual-stress modeling. To establish a point of departure, we briefly 
describe the ADM. 

Let Uj be the fluid velocity, overbars denote linear filtering, and Uj denote the k- times 
successively filtered velocity. The ADM models the exact residual stress 

Rij = uiy,j - UiUj (!) 


as Mij ~ Rij , where 


and 


= ViVj - ViVj 


Vj 


±c k *<r» 


( 2 ) 

(3) 


The quantity v jy termed the deconvolved velocity, approximates uj by defiltering Uj. For LES, 
the coefficients Ck are judiciously chosen so that Vj faithfully restores the low wavenumber 
content of Uj while attenuating high-wavenumber content. Thus, the method is useful for 

‘Research supported in part by NASA Grant NAGl-02033 and by the College of Science and Mathematics, 
James Madison University. 
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LES expressly because it is approximate, not despite its being approximate. The index p 
determines the order of the deconvolution. 

Suppose if (<f) is the transfer function that defines the action of the filter in Fourier space, 
where £ = /cA, k is the wavenumber, and A is the filter width, which parameterizes the filter. 
Provided the filter is invertible (if ^ 0) and |1 - if| < 1, the transfer function has an exact 
power-series inverse, namely 

ir 1 = — -r~r = 1 + (1 - H) + (1 - Hf + ... + {1 - HY + ... (4) 

Truncating the series at finite order p yields the transfer function of the approximate inverse 
if -1 , namely 

if" 1 = 1 + (1 - if) + (1 - if) 2 + ... + (1 - ff) p (5) 

The product if * H~ l is the Fourier-space analog of Eq. 3. By isometry between Fourier and 
physical space, the coefficients C* are determined simply by the binomial theorem (Pascal’s 
triangle); that is, for example, for p = 3, [C 0 > Cj, C 2 , C 3 ] = [4, —6, 4, —1]. 

As Stolz and Adams (1999) note, the ADM (Eq. 2) can be viewed as a generalized 
scale-similarity model (GSSM). It is well known that similarity models are insufficiently dis- 
sipative without secondary regularization (high-order artificial viscosity) and tend to produce 
numerical instabilities. Stabilization is accomplished in the original ADM by the addition of 
a dissipative term to the right-hand side of the momentum equations, namely 

~X(Uj ~ vj) ( 6 ) 

where x is an arbitrary damping parameter, to be determined. Thus, in the ADM, decon- 
volution serves two distinct purposes: 

1. Modeling of the residual stress 

2. Generation of high-order artificial viscosity 

In Fourier space, the operators associated with the two purposes above are, respectively 

1. if* if" 1 

2. l-if*#" 1 

If the filter is based upon a symmetric stencil, the transfer function if is purely real, as are 
H * H- 1 and 1 - if * H~\ 

Let H~ l denote the approximate inverse of order p. Figure 1 presents the transfer function 
of a centered, second-order, parameterized Pade filter, with the nominal cutoff set at £ c = 
7r/2. Also shown are the exact inverse, the approximate inverse for p = 3, and the two 
operators of interest: H * if -1 and 1 — if * H~ l . Strictly speaking, H is not invertible; 
consequently, its exact inverse is unbounded. However, its approximate inverse is well-defined 
for all finite orders p. 

In closing this section, the point to be made is that the “Pascal” coefficients serve well 
in both capacities of the deconvolution: modeling and dissipation. The operator H * H~ l 
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manifests desirable low-pass properties: it is unity at f = 0, is flat in the vicinity of the 
origin (which indicates that low wavenumber content is faithfully recovered by defiltering), 
and drops off rapidly to zero at high wavenumbers. Because symmetric spatial filtering 
induces no phase error, Uj and Vj are aligned in phase, and the operator 1 — H * H ' 1 
is purely real. The product of y(l - H * H~ l ) quantifies the exponential decay rate as a 
function of wavenumber. If 1 — H* H~ l is of one sign (as it is for the Pascal coefficients), and 
X > 0, then secondary regularization is purely dissipative, with the dissipation concentrated 
at high wavenumbers (Fig. 1). 

In the next sections, we describe modifications necessary to adapt approximate decon- 
volution from spatial to time-domain filtering. The resulting temporal variant of the ADM 
will be referred to as the TADM; similarly, LES with with the TADM will be termed TLES 
(temporal LES). 


2 Causal Time-Domain Filtering 


Consider an exponential causal filter expressed in differential form, namely 


dUj \Lj 'ULj 

~dt = A 


(7) 


and parameterized by the temporal filter width A. The transfer function of the continuous 
filter above is 


Hfa A) 


1 

1 4- iuj A 


( 8 ) 


where u> is the dimensional circular frequency and i = i/— I. 

When the filter is implemented by the discrete solution of Eq. 7 with time increments 
At, it is natural to re-parameterize in terms of the filter- width ratio r = A/ At, in which 
case the transfer function becomes 




1 _ 

1 -f irQ, 


(9) 


where 0 = uAt. If the ODE is solved accurately, say, by a high-order numerical method such 
as fourth-order Eunge-Kutta, then the transfer function of the fully discretized problem is 
virtually identical to Eq. 9, so that the continuous transfer function may be used for purposes 
of analysis. 

Whereas spatial filters with symmetric stencils are characterized by purely real transfer 
functions, causal filters have complex transfer functions, because they are by definition biased 
in time, as, for example, in Eq. 9. As a consequence, causal filtering produces both amplitude 
attenuation (desired) and phase error (undesired). The phase shift <f> of the exponential filter 
above is given by 

0(f l;r) = tan -1 — — tan -1 (rQ) (10) 
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Because of the phase error associated with causal filtering and subtle differences be- 
tween the transfer functions of causal and spatial filters, both components of the ADM 
(residual-stress modeling and secondary regularization) require modification when adapted 
for time-domain filtering. This presents two possible approaches to addressing the prob- 
lem: either a single deconvolution for both purposes, or independent deconvolutions, each 
designed for a specific purpose. Here we follow the latter path. Dual deconvolutions re- 
quire some additional storage but no appreciable additional computational effort. The next 
two sections address deconvolution optimizations for residual-stress modeling and secondary 
regularization, respectively. 


3 Residual- Stress Model 

Formally, the residual-stress model of the TADM is identical to Eq 2. Phase error is not 
problematic for the residual-stress model, because both terms of Eq. 2 have the same phase 
relationships. On the other hand, because H (see Eq. 9) is complex for causal filters, other 
problems arise for which adjustments must be made. 

Recall that in Fourier space, the approximate inverse of H is given by 

H~ l = £ C k H k (11) 

k=0 

Insight into the efficacy of the deconvolution can be gained through analysis of the operator 
which is complex. Ideally, for applications to LES, the modulus \H*H~ 1 \ is shaped 
like a spectral (sharp cutoff) low-pass filter. The closer the operator approximates this ideal, 
the better the residual-stress model is likely to be, as implied by Fig. 2, which presents |ifj 
and j H * H~ 1 j relative to the spectral ideal for a deconvolution of order p — 3. The ideal 
exactly restores low? w?avenumber content (resolved scales) while attenuating energy at high 
wavenumbers (unresolved scales) to zero. 

The Pascal coefficients of the ADM are unsuitable for the TADM as shown in Fig. 3. Not 
only is energy at moderate -wavenumber not attenuated, it is amplified. Thus, an optimized 
TADM requires that the coefficients C k be expressly engineered to give H * H~ l “near-ideal” 
properties. 

To illustrate the design process, consider a third-order deconvolution (i.e., p — 3). The 
parameter space consists initially of the four free parameters [Co, C2, C3]. The require- 
ment that the coefficients sum to unity reduces the parameter space by one. That is, 

C 3 = 1.0 - (Co + Ci + C 2 ) (12) 

The remaining coefficients are found by forcing to zero successive derivatives of the modulus 
of H * H~ l , which flattens its graph in the vicinity of the origin. 

There is, however, a further consideration that involves the the zeroth coefficient, Co- 
Unlike spatial filters, such as the Pade filter shown previously, causal filters do not in gen- 
eral vanish identically at the Nyquist frequency (fi = 7r); rather, they tend toward zero 
asmyptotically as Q — > 00, as, for example, in Eq. 9. Because all terms of H~ 1 except the 
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zeroth-order term contain H (see Eq. 11), H~ l — * Co as Cl — oo, and if Co 7^ 0, H * H ~ 1 
tends toward zero as 1/Cl. For example, with Co = 1, H * H~ l — * H as fl — * 00 , as seen 
in Figure 2. However, by setting C 0 = 0, attenuation at high frequencies increases to 1 /Cl 2 , 
which is highly desirable (Fig. 4). From a physical point of view, setting C 0 = 0 implies that 
the deconvolved velocity is re-constructed only from fields that have been low- pass filtered 
at least twice. This ensures that the process of deconvolution is stable at high frequency. 

At this point, only C\ and C 2 remain free parameters. These are determined by setting 
the second and fourth derivatives of the modulus of H * H~ l to zero (as the odd derivatives 
are automatically zero). Algebraic constraint equations can be readily derived by use of 
computer algebra software. In general, the resulting equations are nonlinear and admit 
multiple solutions. Complex solutions, which are discarded because they cannot be readily 
implemented in physical space, occur in conjugate-symmetric pairs. The remaining real 
solutions yield approximate inverses H~ l that share the same moduli but differ in phase. In 
particular, for p — 3, the method described herein yields two real and two comple x solution s, 

from which the following set of real coefficients is selected as optimal: [0.0, \/6, y 4 + 2\/6 — 

2y/6 , 1 — 1/4 4 - 2\/6 + 2 %/ 6 ). Figure 4 presents the moduli of H, Hr 1 , and H * H~ l for 
this coefficient set. The zero-valued derivatives of j H * H~ l \ force its graph to be as flat 
as possible at low frequencies, while the value of zero for Co causes the sharp drop-off at 
higher frequencies. The resulting effect is an amplitude restoration near zero, and a greater 
attenuation at high frequencies (Fig. 4). 

Obviously the method can be extended to higher order deconvolution. 

Experience to date suggests that the temporal residual-stress model is more than ad- 
equate. Deconvolution orders of 2, 3, 4, and 5 have been examined, as have a variety of 
deconvolution coefficients and filter widths. In all cases, the instantaneous and mean resid- 
ual stresses are well defined, with mean components that have qualitatively correct profiles. 
In sum, temporal deconvolution appears to admit a family of robust residual-stress models. 

4 Secondary Regularization 

The purpose of secondary regularization is to impose high-order artificial viscosity at high 
frequencies while leaving low frequencies relatively undisturbed. 

Consider the complex exponential function u(t) = e wt , which satisfies the complex dif- 
ferential equation of a harmonic oscillator, namely 

^ = iuu (13) 

dt 

Secondary regularization functions as a dissipative term for a harmonic oscillator. Accord- 
ingly, consider the model problem 

^ = ium - X (u - H * H~ l u) = [iu + X (H * H - 1 - l)]w 

at 

The additional term, the Fourier-space analog of Eq. 6, imposes exponential damping pro- 
vided 

Re[iw + x{H *H~ l -1)] <0.0 
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which is satisfied iff 

Re\ x (H * H “ a - 1)] < 0.0 
or 

Re\H * H* 1 } < 1.0 (14) 

The exponential envelope (e.g. e Xt ) of the damped harmonic oscillation is governed by the 
decay rate A(O) = x * [He(H * H~ l ) - 1], which in turn is scaled by the damping parameter 

Constraint Eq. 4 is violated whenever the Pascal coefficients are used in secondary regu- 
larization in the TADM. As shown in Fig. 5, the decay rate is positive for some frequencies, 
which leads to unstable exponential growth in time. 

For a given order p, optimal coefficients can be found by setting derivatives of the real 
part of H * H~ l to zero, which, once again, is readily accomplished using computer algebra 
software. For example, for p = 3, there are four coefficients, three of which are free param- 
eters. Optimal coefficients are found by forcing the second, fourth, and sixth derivatives 
of Re[if * H~ 1 } to zero. (Odd-order derivatives are automatically zero.) Specifically, the 
optimal coefficients for third-order temporal secondary regularization are [Co, Ci, C2, C3] = 
[||, ^1], The decay rate for this coefficient set is shown in Fig. 6 (for a x of unity). 

Similarly, an optimal set of coefficients for second-order (p = 2) secondary regularization is 
[C 0) Ci, C 2 ] = [y , j,[]. Figure 6 also presents the decay rate for optimized p = 2 secondary 
deconvolution. 


5 Results 

Figures 7, 8, and 9 compare results of TLES at nominal Re r = 590 to reference channel- 
flow data from DNS by Moser et al (1988). Parameter values for the TLES and refer- 
ence simulations are presented in Table 1. The TLES is conducted with a TADM of or- 
der p — 3 with distinct coefficient optimizations for residual-stress modeling (coefficients 
[0.0,2.4495,-1.9159,0.4664]) and secondary regularization (coefficients [f§, ^ , f, ^]), as 
discussed in the previous two sections. The flow is initialized to a laminar state, randomly 
perturbed. Mean quantities are computed by averaging over 500 < t < 1500. Spatial reso- 
lution for the TLES is commensurate with that of the spatial LES of channel flow with the 
ADM by Stolz et al (1999a), whose results are also presented in Table 1. 


Case 

N x x N y x N z 

At 

-^ e hulk 

RCrJ- 

model 

(d x) 

KMM 

384 x 384 x 257 

— 

10935 

587 

none 

NA 

SAK 

48 x 64 x 65 

— 

10935 

574 

ADM 

NA 

TLES 

48 x 64 x 65 

0.04 

10935 

595 

TADM 

(8,1.0) 


Table 1: Parameter values of channel-flow TLES at nominal Re T — 590 and reference simu- 
lations of Moser et al (1988) (KMM) and Stolz et al (1999a) (SAK). 
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The present results of TLES with the TADM are highly encouraging. With optimized 
secondary regularization, the TADM is sufficiently robust to survive the enormous spike in 
turbulent kinetic energy k during transition (Fig. 7). The TADM adjusts accordingly during 
transition as indicated by the spike in k R (the trace of the modeled residual-stress tensor 
My) during the time interval 50 < t < 150. The computed value of Re r = 595 deviates only 
1.3 percent from the reference value of 587, an overshoot commensurate with the undershoot 
of the results of Stolz et al (1999a) using the ADM. Meanflow profiles axe quite good near 
the wall but deviate somewhat from the DNS results near center channel (Fig. 8). Reynolds- 
stress profiles agree well qualitatively and quantitatively with DNS reference profiles (Fig. 9), 
with significant contribution to the Reynolds stress from the mean residual stress. 

Table 2 summarizes results from several simulations of the present method (TLES) at 
nominal Rer = 180, for which, for brevity, no figures are presented. Results from a reference 
LES of Stolz et al (1999a) using the spatial ADM are also provided. All simulations were 
conducted at a spatial grid resolution of 32 x 32 x 33, whereas the reference DNS of Moser 
et al (1988), for which Re T = 178.1, was conducted at resolution 128 x 128 x 129. 

The results of TLES (with optimized coefficients) at nominal Re T = 180 are somewhat 
disappointing relative to the results at nominal Re r — 590. All trial cases represent im- 
provement over the “no-model” case (TLESl80a), in which both the residual-stress model 
and secondary regularization are turned off. However, computed Re r is nearly 10 percent in 
error for Case TLESl80b, which has the same deconvolution coefficients as for the Re r = 590 
simulation presented previously. It should be noted that, although Re T = 590 represents a 
more computationally intensive calculation than Re T = 180, Re r = 180 may represent the 
more severe test of the TADM, because barely turbulent flow is in general more anisotropic 
than highly turbulent flow. Comparison of Cases TLBSl80b-c reveals the solution to be 
relatively insensitive to the damping coefficient y, although slightly improved by higher dis- 
sipation. Because the damping coefficient is arbitrary, this is a positive result; it would be 
highly undesirable for the solution to depend sensitively on an arbitrary value. (Note that 
the ADM also involves an arbitrary coefficient for secondary regularization, to which the 
solution is also relatively insensitive.) 

Case TLESl80d results from mixed deconvolution coefficients, where optimized third- 
order ( p — 3) coefficients are exploited for the residual-stress model, but second-order (p = 2) 
optimization is used for secondary regularization. Note that both deconvolutions have the 
same number of free parameters (namely two), which in each case is used to force appro- 
priate second and fourth derivatives to zero at the origin. That Case TLESl80d represents 
further improvement suggests that higher order deconvolution is not necessarily better with 
regard to artificial dissipation. Here we interject some speculation. So called mixed models 
of residual stress blend a scale-similarity model with a Smagorinsky-like dissipative term 
(secondary regularization) as an attempt to fix the inadequate dissipation of a stand-alone 
scale-similarity model. The result is typically overly dissipative, in our view, because the 
dissipation comes in at second order. Our current thinking is that the dissipative term needs 
to be of higher order than the second derivative (at which the physical viscosity is active), 
but not necessarily of very high order. 

In conclusion, the present addendum outlines the constraints that should be imposed 
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Case 

c 2 

X 

Re T 

model 

TLES 180a 

N/A 

N/A 

203.0 

none 

TLES 180b 

[2.1875, -1.8125, 0.7500, -0.1250] 

1.0 

195.4 

TADM 

TLES 180c 

[2.1875, -1.8125, 0.7500, -0.1250] 

2.0 

193.5 

TADM 

TLESl80d 

[1.8750, -1.1250, 0.2500, 0.0] 

1.0 

189.8 

TADM 

SAKl80a 

N/A 

N/A 

173. 

ADM 


Table 2: Results for various trials at nominal Re r = 180, with order p = 3, filter- 
width ratio r — 8, ^«b\alk = 2800, vector of coefficients in residual-stress model C x = 
[0.0, 2.4495, —1.9159,0.4664], and vector of coefficients in secondary regularization Ca as in- 
dicated. 

to optimize temporal deconvolution in the TADM for the dual purposes of residual-stress 
modeling and secondary regularization. We have focused on third-order (p = 3) deconvo- 
lutions, but the optimization approach is applicable to any order p. Furthermore, we have 
experimented very little with the filter-width ratio r, which remains an essential parameter 
of the TADM. Consequently, present coefficients should not be considered ideal in any sense. 
In sum, TLES with the TADM appears to be a promising alternative to LES that affords 
considerable opportunity for refinement. 
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Figure 2: Transfer functions for H (dashed), H 1 (dotted), H * H 1 (dashed and 

dotted), and spectral function (solid) for a third-order deconvolution with coefficients 
[1.0, 0.559, -0.784, 0.225] and r = 2. 



CU 71 


Figure 3: Transfer functions for H (dashed), H~* (dotted), H*H~ l (dashed and dotted), and 
spectral function (solid) for a third-order deconvolution with Pascal coefficients [4, —6, 4, — 1] 
and r = 2. 
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Figure 4: Transfer functions for H (dashed), H 1 (dotted), H * H 1 (dashed and dot- 
ted), and a spectral function (solid) for third-order deconvolution with optimized coefficients 
[0.0,2.4495,-1.9159,0.4664] and r = 2. 



Figure 5: Decay rate (A) for third-order deconvolution with Pascal coefficients [4, —6, 4, -1], 
r = 2, and x — 1-0. 
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Figure 6: Decay rates (A) for third-order and second-order deconvolutions with optimized 
coefficients [2.1875, —1.8125, 0.75, —0.125] and [1.875,-1.125,0.750], respectively, r — 2, and 
X = 1.0. LEGEND: p = 3 (solid); p = 2 (dashed). 



Figure 7: Evolution of Re T and k for Case TLES590. LEGEND: solid line (instantaneous); 
dashed line (mean); dashed-dotted line (k); dotted line (/cr). 
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Figure 8: Mean streamwise velocity profile < U\ > vs. z and its wall-normal derivative for 
Case TLES590 relative to reference DNS data of Moser et al. LEGEND: solid line (KMM); 
dashed line (TLES590); symbols (TLES590, wall units). 
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